function LEs = LEs(state0, mu, k, nSteps)
    % 计算 MCLM 模型的 Lyapunov 指数
    % LE1: 代表最大的李指数，反映了系统对初始扰动中最剧烈的发散方向，通常决定了系统是否具有混沌行为
    %      若 LE1 为正，则系统呈现混沌或超混沌状态; 若 LE1 为零或负，则可能处于周期性或稳定状态
    % LE2: 系统在与 LE1 正交的方向上扰动的平均扩展或收缩速率
    % LE3：最小的李指数，它表示系统在剩余方向上的扰动行为
    % LE1 + LE2 + LE3 的总和反映了系统体积的收缩或扩展情况。如果总和为负，则系统具有吸引性质
    
    d = 3;              % 系统维数
    Q = eye(d);         % 初始化切向量（切空间基向量）
    LE_accum = zeros(1, d);  % 用于累计各方向的对数伸缩因子
    
    state = state0;     % 当前状态初始化
    
    for n = 1:nSteps
        % 更新状态（调用 mclm 函数）
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
        
        % 计算当前状态的 Jacobian 矩阵
        J = jacobian(state, mu, k);
        
        % 更新切空间基向量并正交化
        Q = J * Q;
        [Q, R] = qr(Q);
        
        Q = Q * diag(sign(diag(R))); % 避免方向翻转

       % 累加 R 矩阵对角元素的绝对值对数
        LE_accum = LE_accum + log(abs(diag(R)))';
        
    end
    
    LEs = LE_accum / nSteps;  % 计算平均 Lyapunov 指数
end
